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Numerical  simulations  were  performed  with  a  0.2 

_ d*gpe^r  two  layer,  free  surface,  primitive  equation,  ocean 

circulation  model  of  the  Gulf  of  Mexico  with  a  realistic 
bottom  topography  and  coastline.  As  in  previous 

simulations,  with  idealized  bottom  topography  and  coastline, 
realistic  Loop  Current  behavior  could  be  obtained  without 
wind  forcing  and  with  constant  inflow  through  the  Yucatan 
Straits.  However*  simulations  with  both  wind  and  port 
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forcing  showed  significantly  more  variation  than  those  with 
port  forcing  alone. 

A  similar  two  layer  model  f  has — -been  set  up  for  a 
North  Atlantic  region, ^  from  20N,  82W  to  48N,  30W,  on  a  0.25 
by  0.2  deg'ree^grid.  4  '?;*  a  P  •>  >  . 

-A- J'BlarO-k—  frcr**  -  mult  igr  id  package  fras~b-e^n/  prepared 
for  solving  finite  difference  Helmholtz's  equations  with 
staggered  grid  Neumann  boundary  conditions  in  nonrectangular 
regions.  In  semi -implicit  free  surface  primitive  equation 
ocean  models,  it  is  competitive  with  other  powerful 
iterative  methods  but  slower  than  the  direct  papacitance 
Matrix  Technique.  l( , . ,  j 

Several  wind  stress  data  sets  *~ha-ve  been'sprepared 
for  use  in  ocean  models,  based  on  global  analysis  products 
from  FNOC  and  NMC,  regional  analysis  products  from  FNOC,  and 
ship  observations.  These  data  sets  have  been  used  to  drive 
ocean  circulation  models  of  the  World  Ocean,  the  Indian 
Ocean,  the  North  Atlantic,  and  the  Western  Mediterranean 
Sea.  Y ■ 

The  capabilities  of  a  statistical  model  to  nowcast 
subsurface  pressure  information  from  free  surface  data 
obtained  by  a  hypothetically  perfectly  accurate  altimeter 
was  explored.  Pointwise  linear  regression  of  upper  to  lower 
layer  pressure  did  not  lead  to  successful  lower  layer 
predictions.  Pointwise  multilinear  regression  was  more 
successful,  particularly  in  predominantly  barotropic 
regions.  However  it  did  not  do  well  in  regions  of  very 
strong  baroclinic  instability.  Statistical  prediction 
appears  to  be  a  promising  technique  for  initializing 
dynamical  ocean  circulation  prediction  models. 
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GULF  OF  MEXICO  MODELING 


This  is  a  continuation  of  modeling  in  the  Gulf  of 
Mexico  previously  performed  by  NORDA  (Hurlburt  and  Thompson, 
1980  and  1982).  The  same  two  layer,  primitive  equation, 
free  surface,  hydrodynamic  model  with  topography  confined  to 
the  lowest  layer  is  used  as  before,  but  with  the  addition  of 
'realistic'  (rather  than  idealized)  coastline  geometry  and 
bottom  topography.  The  ocean  model  contains  many  innovative 
features  and  the  reader  is  referred  to  its  detailed 
discussion  in  Hurlburt  and  Thompson  (1980).  In  particular, 
Section  2  (pp  1613~ 161M)  gives  the  model  equations  and 
Appendix  B  (pp  1647-1650)  describes  the  numerical  design  of 
the  model.  Since  that  time  the  capability  to  handle  general 
basin  geometry  has  been  added  but  this  does  not  affect  the 
description  in  any  major  way.  Wind  forcing  is  treated 
identically  to  interfacial  and  bottom  stress  terms,  i.e., 
wind  stress  appears  directly  as  an  additive  term  in  the 
momentum  equation  [see  p  1614  of  Hurlburt  and  Thompson 
(1980)]. 

In  terms  of  'realism'  Experiment  9  was  the  most 
successful  port  driven  Gulf  of  Mexico  numerical  simulation 
produced  under  this  contract.  The  model  was  driven  from 
rest  to  statistical  equilibrium  solely  by  a  steady  inflow 
through  the  Yucatan  Straits  which  was  compensated  by  outflow 
through  the  Florida  Straits.  The  model  parameters  were: 

•  upper/lower  layer  inflow  transport  -  26/4  Sv, 

•  wind  stress  -  0, 

•  horizontal  eddy  viscosity,  A  -  300  m2/sec, 

•  grid  spacing  -  25  by  25  km, 

•  reference  layer  thicknesses,  HI  -  200m  and 

H2-3400m, 

•  minimum  depth  of  bottom  topography  -  500m, 

•  beta,  df/dy  -  2 *  1 0 ~ 1 1  m~*  see-1, 

•  Coriolis  parameter  at  the  southern  boundary, 
f  -  5  x  1 0" *  sec-1 , 


•  gravitational  acceleration,  g  -  9.8  m/sec2, 

•  reduced  gravity,  g'  -  . 03  ( HI +  H2) /H2  m/sec2, 

•  interfacial  stress  -  0, 

•  coefficient  of  quadratic  bottom  stress  -  .002,  and 

•  time  step  -  1  hour. 

Figure  1  compares  'instantaneous'  upper  ocean  flow 
patterns  just  before  an  eddy  is  shed  by  the  Loop  current  (a) 
from  the  numerical  model  and  (b)  from  observations  by 
Leipper  (1970).  The  ability  of  the  model  to  simulate 
observed  features  is  clearly  demonstrated  by  this 
comparison,  which  is  remarkable  given  the  simplicity  of  the 
model  forcing.  However  some  discrepancies  remain,  for 
example  the  eddy  has  not  penetrated  a3  far  into  the  Gulf  and 
is  more  intense  than  that  shown  in  the  observations.  Waves 
can  be  seen  moving  around  the  wall  of  the  Loop  Current  in 
both  the  model  and  the  observations,  but  in  the  model  they 
are  at  the  limit  of  resolution  and  therefore  unrealistically 
large.  Moreover  in  the  Gulf  the  waves  are  more  pronounced 
on  the  eastern  wall  of  the  Loop  and  can  form  strong  cold 
intrusions  that  may  contribute  to  the  eddy  shedding  process 
(Vukovitch  and  Maul,  1984).  As  shed  eddies  propagate 
westward  (Figure  2a)  the  model  spontaneously  develops  a 
counter  rotating  vortex  pair  (Figure  2b),  a  structure 
repeatedly  observed  in  the  Western  Gulf  (Figure  3).  The 
roles  of  the  wind  and  the  Loop  current  eddies  in  the 
formation  of  this  structure  have  been  a  matter  of  some 
controversy  (Merrell  and  Morrison,  1981).  Although  wind 
forcing  was  not  present  in  this  simulation  a  major  role  for 
winds  has  not  yet  been  ruled  out.  After  spin-up  the 
experiment  sheds  an  eddy  once  every  390  days  and  the  eddy 
shedding  cycles  are  very  similar. 

Experiment  28  is  identical  to  Experiment  9,  except 
for  the  grid  spacing  (20  by  22,  rather  than  25  by  25  Km)  and 
the  Coriolis  parameter  (4.5x10“*,  rather  than  5x10“*  see-1). 
It  gives  very  similar  results.  Experiment  31  uses  the  same 


model  parameters  as  Experiment  28,  but  has  no  port  forcing 
being  driven  instead  by  seasonal  climatological  winds  based 
on  ship  observations  (Elliot,  1879).  Linear  interpolation 
in  time  was  used  between  the  seasonal  fields  to  produce  the 
wind  stress  at  each  time  step.  The  experiment  reaches  an 
annual  steady  state  cycle,  winds  and  associated  circulations 
for  each  season  are  shown  in  Figures  4  to  7  (note  the 
contour  interval  is  12.5  m  rather  than  the  more  usual  25  m) . 
Upwelling  in  the  south  west  Gulf  in  the  fall  and  winter,  and 
a  persistent  anticyclonic  'Texas*  coastal  current  are  the 
major  features  of  the  simulation. 

Experiment  34  combines  the  wind  and  port  forcings 
of  Experiments  28  and  31.  Model  parameters  were: 

•  upper/lower  layer  inflow  transport  -  26/4  Sv, 

•  wind  stress  from  seasonal  climatology  based  on  ship 
observations, 

•  horizontal  eddy  viscosity,  A  -  300  m2/3ec, 

•  grid  spacing  -  20  by  22  km  (0.2  by  0.2  degrees), 

•  reference  layer  thicknesses,  HI  •  200m  '  d 

H2-3300m, 

•  minimum  depth  of  bottom  topography  -  500m, 

•  beta,  df/dy  -  2 *  1 0 “ 1 1  m-1  sec-1, 

•  Coriolis  parameter  at  the  southern  boundary, 
f  -  4.5  x  10-5  sec-1, 

•  gravitational  acceleration,  g  -  9.8  m/sec*, 

•  reduced  gravity,  g’  -  .03 ( HI +H2) /H2  m/sec2, 

•  interfacial  stress  -  0, 

•  coefficient  of  quadratic  bottom  stress  *  .002,  and 

•  time  step  -  1.5  hours. 


Seasonal  climatological  winds  obviously  cannot 
represent  most  of  the  wind  variability  and  persistence. 
Even  so  the  addition  of  wind  forcing  increases  the 
variability  of  the  Loop  Current  system,  including  the  eddy 
shedding  period  and  eddy  path.  For  example  Figure  8 
compares  360  model  days  from  Experiments  28  and  34  (which 


are  identical  except  that  Experiment  28  has  no  wind 
forcing).  From  these  snapshots,  taken  every  90  days,  there 
is  little  difference  between  the  two  experiments.  But  if 
Experiment  34  is  sampled  every  20  days,  as  in  Figure  9,  it 
is  apparent  that  two  eddies  were  shed  in  the  space  of  about 
one  year.  Figure  9  also  shows  that  the  circulation  pattern 
in  the  Western  Gulf  can  change  very  rapidly  at  times. 
Figure  10  shows  the  mean  interface  deviation  and  its 
variability  for  experiments  with  wind  forcing  only,  port 
forcing  only  and  with  wind  plus  port  forcing.  This 
demonstrates  that  even  in  the  mean  the  interaction  of  wind 
and  port  forcing  is  not  linear,  i.e.,  the  mean  of  the  dual 
forcing  experiment  is  not  the  sum  of  the  other  two  means. 
The  variability  is  increased  in  the  dual  forcing  case, 
particularly  in  the  central  and  western  Gulf. 

More  recent  work  has  largely  superseded  the 
experiments  described  above.  In  particular  a  ratio  of 
20/10  Sv  for  upper/lower  layer  transport  is  now  preferred 
over  the  26/4  Sv  used  here,  and  wind  driven  experiments  have 
been  performed  with  wind  data  based  on  the  Navy's  twelve 
hourly  surface  pressure  analysis  from  1967  to  1982 
(Wallcraft,  1985).  The  differences  between  the  seasonal 
means  from  the  12  hourly  data  and  the  ship  observation  data 
is  significant,  particularly  in  the  south  west  Gulf  (Rhodes, 
Thompson  and  Wallcraft,  1985).  Therefore  emphasis  should 
not  be  placed  on  the  details  of  the  wind  driven  experiments 
described  above. 


II 


NORTH  ATLANTIC  MODELING 


The  aim  of  this  modeling  effort  was  to  set  up  a 
finite  depth  ocean  model  of  a  large  ocean  basin.  Detailed 
modeling  experiments  were  not  performed.  The  major 
difference  between  the  small  basin  models  (of  the  Gulf  of 
Mexico,  for  example)  and  a  North  Atlantic  model  is  that  the 
beta-plane  approximation  no  longer  applies  so  spherical 
coordinates  must  be  used.  Explicit  reduced  gravity 
versions  of  large  basin  models  already  existed  (NORDA's 
World  Ocean  Model,  for  example),  and  the  major  difficulty  in 
producing  the  finite  depth  model  was  the  development  of  the 
correct  Helmholtz  solver  for  the  spherical  semi-implicit 
equations.  At  the  same  time  the  two-layer  finite  depth 
model  was  generalized  to  n-layers,  for  both  the  beta-plane 
and  for  spherical  coordinates.  Thi3  involved  the 
replacement  of  precalculated  constants  for  transformation 
into  modal  space,  with  numerically  derived  values. 

The  model  region  is  from  20N,  82W  to  48N,  30W  and 
all  boundaries  are  closed  with  free  slip  conditions  being 
applied  to  the  east  and  south  boundaries.  The  grid  spacing 
is  0.25  by  0.2  degrees,  i.e.,  delta-x  varies  from  26  km  (at 
20N)  to  19  km  (at  48N)  and  delta-y  is  about  22  km.  This  is 
the  smallest  closed  region  that  is  viable  for  studying  the 
Gulf  Stream.  A  larger  region  (and  finer  grid)  would  be 
preferable,  but  would  require  a  faster  machine  with  more 
main  memory  than  the  CRAY  X-MP  1200  used  for  these 
experiments. 

Figure  11  shows  the  coastline  and  bottom  topography 
used  by  the  model.  The  topography  is  from  SYNBAPS,  the 
minimum  depth  is  1000  m  and  it  has  been  smoothed  twice  by  a 
nine  point  real  smoother.  Figure  12  is  a  plot  of  interface 
deviation  from  the  seventh  year  of  an  experiment  forced, 
from  rest,  with  monthly  climatological  winds.  A  meandering 
Gulf  Stream  has  been  spun  up  by  this  simulation,  however  its 


separation  point  from  the  coast  and  penetration  distance  are 
not  ’realistic'.  The  simulation  was  only  performed  to  check 
out  the  ocean  model,  it  is  expected  that  a  full  scale 
modeling  effort  in  the  North  Atlantic  will  be  started  by 
NORDA  in  the  near  future  to  take  this  work  further. 


III. 


MULTIGRID  SOLVER. 


Multigrid  techniques  are  a  relatively  new  set  of 
methods  for  greatly  accelerating  the  convergence  of  the 
classical  iterative  techniques  for  the  solution  of  elliptic 
partial  differential  equations  (Brandt,  1982).  Multigrid 
iteration  involves  a  sequence  of  grids,  each  with  about  four 
times  as  many  nodes,  or  computational  points,  as  its  coarser 
predecessor.  On  each  grid  (except  the  coarsest),  high 
frequency  error  components  are  reduced  by  a  classical 
iterative  method  and  low  frequency  errors  by  extrapolation 
from  the  solution  of  a  related  problem  on  the  next  coarsest 
grid.  The  methods  great  speed  is  based  on  the  fact  that 
many  iterative  procedures  reduce  high  frequency  errors  very 
rapidly,  and  also  on  the  relative  cheapness  of  solutions  on 
coarse  grids. 

In  NORDA's  layered  semi-implicit  primitive  equation 
free  surface  ocean  circulation  models  the  solution  of 
several  finite  difference  Helmholtz 's  equations  are  required 
to  relatively  high  accuracy  (about  5  or  6  places  of 
decimals)  each  time-step  (one  per  layer).  These  equations 
have  staggered  grid  Neumann  boundary  conditions,  and  are 
over  the  general  irregular  regions  formed  by  ocean  basins 
(not  over  a  rectangle,  which  would  be  computationally 
simpler).  Classical  multigrid  methods  could  be  applied  to 
such  problems,  but  the  irregular  geometry  and  Neumann 
boundary  conditions  give  rise  to  bookkeeping  and  algorithmic 
difficulties.  An  alternative  is  to  use  a  "black  box" 
multigrid  method  (Dendy,  1982),  which  acts  on  the  general 
finite  difference  form  of  the  equations  and  makes  minimal 
assumptions  about  the  underlying  continuum  problem. 

A  public  domain  black  box  multigrid  solver  (Dendy, 
1982)  has  been  modified  for  NORDA's  ocean  models.  The 
original  code  was  for  general  coefficient  finite  difference 
9-point  elliptic  equations  in  a  rectangle  with  discontinuous 


coefficients.  This  code  applies  equally  to  five-point 
Helmholtz's  equation  with  Neumann  boundary  conditions  over  a 
nonrectangular  region,  since  the  irregular  boundary  is 
equivalent  to  a  discontinuity.  Black  box  multigrid  always 
uses  a  nine-point  operator  on  the  coarser  grids,  even  when 
the  finest  operator  is  five-point.  On  the  other  hand  the 
irregular  boundary  needs  no  special  coding  on  the  coarse 
grids.  Modifications  to  the  original  code  include: 

a)  Uses  red-black  (four  color  for  nine-point 

equations)  pointwise  Gauss  Seidel  as  the  underlying 
iterative  scheme,  because  it  is  highly 

vectorizable . 

b)  Modified  finest  grid  code  to  handle  the  required 

Helmholtz's  equations  only.  This  saves  both 

storage  and  time  over  the  general  version. 

c)  Fully  vectorized  the  solution  phase  of  the  solver, 
the  preprocessing  stage  has  not  been  vectorized 
since  it  is  only  called  at  the  beginning  of  a  model 
run. 

Total  storage  requirements  are  3.3  meshes  per 
equation  for  constants,  and  2  workspace  arrays.  Table  1 
gives  times,  on  a  CRAY  X-MP  1200,  for  the  multigrid  code  and 
for  a  fast  direct  rectangular  Helmholtz  solver  (based  on  a 
FFT)  over  various  sized  rectangles.  The  multigrid  code  is 
equally  as  fast  on  nonrectangular  regions.  The  FFT  based 
solver  can  be  extended  to  general  regions  via  the 
Capacitance  Matrix  Technique  (Buzbee,  et.  al.,  1980)  but 
solve  times  are  two  to  three  times  as  long  and  a  large 
amount  of  storage  can  be  required. 


TABLE  i  -  GRID  SIZE  vs.  SOLVE  TIME  ana  TIME  VLR  NODE, 

FOR  TWO  HELMHuLxZ  SOLVERS  oN  A  CHAT  X  -  ME  ISO 
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Notes  : 

1)  All  multigrid  times  are  averages  to  reduce  the 
original  relative  error  by  10~*.  This  is  probably 
a  worse  case  for  ocean  models,  where  typical 
required  error  reduction  is  about  10~*. 

2)  Multigrid-1  is  a  case  with  a  small  Helmholtz 

parameter,  and  raultigrid-2  a  case  with  a  large 
parameter.  These  correspond  to  external  and 

internal  gravity  modes  respectively  in  the  ocean 
model. 

3)  Multigrid-1  takes  about  5  to  7  major  iterations, 
Multigrid-2  about  2  to  3  iterations  for  10-# 
improvement  in  accuracy. 

4)  FACR(1)  is  a  fast  direct  solver  that  only  works  for 

rectangles,  minimum  direct  solve  time  over 

irregular  regions  is  2  to  3  times  longer. 


IV 


BOTTOM  TOPOGRAPHIES 


The  global  Synthetic  Bathymetric  Profiling  System 
( SYNBAPS)  data  set  is  uniquely  suited  to  the  generation  of 
ocean  model  bottom  topographies.  It  differs  from  the 
earlier  Northern  Hemisphere  version  in  covering  99  percent 
of  the  world  ocean,  and  in  being  organized  on  a  geographic 
grid  (five  minute  grid  spacing)  rather  than  on  a  meridional 
parts  grid  (Van  Wyckhouse,  1973). 

A  VAX/VMS  specific  FORTRAN  77  package  has  been 
developed  that  extracts  bottom  topography  for  NORDA's  ocean 
models  from  the  SYNBAPS  data  set.  The  capability  to  set  the 
coastline  at  any  desired  depth  below  sea  level,  and  to  cut 
off  the  continental  shelf  at  any  depth  is  included.  The 
coastline  capability  is  required  because  setting  the  coast 
at,  say,  100  m  below  sea  level  might  be  appropriate  in  some 
cases  (particularly  for  reduced  gravity  models).  The 
continental  shelf  cut  off  is  required  because  bottom 
topography  must  be  constrained  to  the  lowest  layer  in 
layered  models,  typically  the  shallowest  topography  might  be 
at  500  m  or  more  below  sea  level.  SYNBAPS  can  give  rise  to 
unacceptable  coastlines,  in  isolated  areas,  either  because 
SYNBAPS  has  no  data  there  or  because  of  interpolation 
difficulties.  Such  areas  must  be  cleaned  up  by  hand,  and  as 
an  aid  to  their  isolation  and  correction,  graphics  display 
software  and  interactive  node  modification  procedures  have 
been  included  in  the  package.  Nine  point  real  smoothers  are 
also  provided,  since  two  delta-x  waves  are  undesirable  in 
model  bottom  topographies. 

Figure  13  shows  the  0.1  degree  Gulf  of  Mexico 
topography,  as  originally  generated  by  the  package  from  the 
raw  SYNBAPS  data.  To  reach  this  stage  requires  about 
fifteen  minutes  (wall  clock  time),  from  entering  the  package 
to  obtaining  a  hardcopy  of  the  plot.  The  coastline  is  taken 
as  the  10  m  depth  contour,  and  all  topography  shallower  than 


500  m  is  set  to  500  m  to  ensure  that  the  topography  is 
confined  to  the  lowest  layer  of  the  ocean  model.  Note  that 
there  is  no  SYNBAPS  data  available  in  the  rectangle  bounded 
by  (95W,  29N)  and  (90W,  30N),  there  are  several  regions  in 
the  world  ocean  where  SYNBAPS  has  no  data  for  small  areas 
near  the  coast.  This  is  apparently  to  reduce  the  size  of 
the  data  base,  since  these  regions  can  only  be  included  if 
large  areas  over  land  are  also  present. 

Figure  14  shows  the  0.1  degree  Gulf  of  Mexico 
topography  modified  for  use  in  the  ocean  model.  The 
coastline  has  been  extensively  edited,  primarily  to  improve 
to  look  of  graphical  output  rather  than  for  modeling 
constraints  (the  actual  model  coastline  is  largely  that 
shown  in  the  previous  figure).  The  Caribbean  Sea  is  not 
part  of  the  model  region,  and  has  therefore  been  set  to  land 
(although  the  Caribbean  coastline  is  still  shown  on  plots). 
Finally  a  nine  point  real  smoother  has  been  applied  to  the 
topography  twice  to  remove  small  scale  variations. 
Production  of  this  version  took  about  one  day,  it  is  of 
about  average  size  ( 1 80  by  130)  but  has  far  fewer  problems 
with  islands  and  narrows  than  some  other  regions. 

The  SYNBAPS  extraction  package  can  only  be  accessed 
interactively,  via  a  DCL  command  file.  It  provides 
sufficient  information  on  the  screen  to  allow  anyone 
familiar  with  the  target  ocean  model  code  to  produce  a 
topography  field  for  a  new  region  without  difficulty. 
Topographies  have  been  created  for  ocean  models  of:  the 
World  Ocean,  the  North  Atlantic  Ocean  (Figure  11),  the 
Indian  Ocean,  the  Gulf  of  Mexico,  and  the  Mediterranean  Sea. 
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V 


WIND  STRESS  DATA  SETS 


The  following  wind  data  sets  have  been  used  to 
compute  wind  stress  fields,  which  in  turn  have  been  or  will 
be  implemented  in  the  World  Ocean  Model  or  other  single¬ 
basin  ocean  models  currently  being  developed  by  NORDA  Code 
323.  In  all  cases,  the  following  equations  were  used  to 
compute  the  x-  and  y-components  of  the  wind  stress,  t: 

tx  -  PCD|tf|u,  xy  -  pCD|tf|v, 

where  p  is  the  air  density  in  g/cm*,  CD  is  the  drag 
coefficient,  for  which  the  value  0.002  was  used,  | V j  is  the 
magnitude  of  the  wind  velocity,  u  is  the  x-component  of  the 
wind  velocity,  and  v  is  the  y-component  of  the  wind 
velocity.  The  exception  was  the  Hellerman  data,  which  was 
received  in  the  form  of  u  and  v  wind  stress  components.  The 
spatial  interpolation  scheme  utilized  was  a  repeated 
application  of  a  one-dimensional  Bessel  interpolation  of  the 
third  order  involving  the  sixteen  values  closest  to  the 
desired  grid  point. 

NOGAPS 


Two  products  of  the  Fleet  Numerical  Oceanography 
Center  (FNOC)  Naval  Operational  Global  Atmospheric 
Prediction  System  (NOGAPS)  analysis  are  the  planetary 
boundary  layer  u-  and  v-components  of  the  wind  (cm/s)  on  a 
1 44  by  73,  2.5°  by  2.5°  grid  that  extends  from  90S  to  90N 
and  6QE  to  57. 5E.  These  wind  components,  available  at  0000Z 
and  1200Z  synoptic  times,  for  the  year  1983  (the  first 
complete  year  available)  were  retrieved  from  FNOC’s  archives 
for  use  in  developing  the  mean  u-  and  v-components  of  the 
wind  stress  for  each  month  of  that  year.  The  procedure  was 
as  follows.  First,  the  wind  stress  components  (dynes/cm2) 


at  each  grid  point  were  computed.  Then  the  monthly  averages 

of  these  components  were  calculated.  Before  storing  the 

results  on  tape  at  FNOC,  the  entire  grid  was  shifted 

westward  so  that  180W  formed  the  western  edge  and  an 

additional  column  was  added  to  the  eastern  edge  of  the  grid, 

a  duplicate  of  the  column  on  the  western  edge.  Thus,  the 

grid  for  t  and  x  now  extends  from  180W  to  180W  in  the  x- 
x  y 

direction. 

The  monthly-mean  wind  stress  components,  derived 
from  the  NOGAPS  u-  and  v-components  of  the  wind,  were  stored 
on  tape  at  FNOC  in  order  to  facilitate  implementation  of 
these  fields  in  the  World  Ocean  Model,  which  is  run  on  the 
FNOC  CYBER  205  vector-processing  computer.  One  additional 
step  before  utilization  by  the  model  was  transfer  of  the  xx 
and  Ty  monthly-mean  fields  to  the  CYBER  205  after  spatial 
interpolation  to  the  243  by  111,  1.5°  by  1.25°  World  Ocean 

Model  grid.  The  x-  and  y-components  of  the  wind  stress  were 
staggered  such  that  the  southwest  corners  of  the  xx  and  xy 
grids  were  at  66.25S, 179.2W  and  65. 625S, 1 71 . 25W, 
respectively.  (There  was  some  overlap  in  the  x-direct ion) . 
Figure  15  shows  the  monthly-mean  wind  stress  curl  for 
January,  April,  July,  and  October  of  1983. 

Hellerman  Wind  Stress  Climatology 

This  data  set  consists  of  the  climatological 

monthly-mean  wind  stress  components,  x  and  x  ,  over  the 

x  y 

world  ocean  on  a  180  by  90,  2°  by  2°  grid  that  extends  from 

89S  to  89N  and  IE  to  359E.  Hellerman  (1967)  computed  these 
wind  stress  components  using  wind  rose  data  and  the  above 
equations,  where 

Cd-2.43  x  1 0“3  for  windspeeds  of  at  least  Beaufort  force  7 
and 
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(  .00220  +  1.136)  x  10“’ 
1.18  x  10~‘ 
(-.0028#  +  1.124)  x  10~J 


0  >  20° 

-20°  <  0  <  20° 
0  <  20° 


These  data  were  interpolated  and  shifted  westward 
onto  the  145  by  73.  2.5°  by  2.5°  grid  used  for  the  NOGAPS 
1983  wind  stress  components  to  facilitate  generation  of 
graphics  using  an  existing  plot  routine  and  comparison  of 
the  two  data  sets.  The  data  were  then  interpolated  to  the 
243  by  111,  1.5°  by  1.25°  World  Ocean  Model  staggered  u-  a^d 
v-grids  and  stored  on  disk  on  the  FNOC  CYBER  205  for  use  in 
the  World  Ocean  Model. 

The  x  and  t  wind  stress  components  were 
x  y 

interpolated  to  the  Indian  Ocean  and  North  Atlantic  Model 
staggered  grids  in  a  similar  manner.  The  Indian  Ocean  Model 
area  extends  from  10S  to  27N  and  35E  to  106E  and  the 
southwest  corner  of  the  u-component  grid  is  at  10S,  38. IE 
and  that  for  the  v-component  grid  is  at  10. IS,  38E.  An 
enlarged  version  of  the  model  extends  from  25S  to  27N  and 
38E  to  106E  so  that  the  lower  left  point  on  the  u-component 
grid  is  at  65S,  38. IE  and  the  same  point  on  the  v-component 
grid  is  at  65.1 S,  38E.  For  both  model  regions,  the  u-  and 
v-wind  stress  components  have  been  stored  on  a  0.2°  by  0.2° 
grid  and  a  0.4°  by  0.4°  grid.  In  the  case  of  the  North 
Atlantic  Ocean,  the  model  region  extends  from  20N  to  48N  and 
82W  to  30W  and  the  wind  stresses  have  been  stored  on  a  20 9 
by  136,  .25°  by  .2°  grid.  The  u-  and  v-grids  are  staggered 
such  that  the  lower  left  point  of  the  tx  wind  stress 
component  is  at  20N,  81 .875W  and  that  of  the  xy  wind  stress 
component  is  at  20. IN,  82W.  Figure  16  is  the  same  as  Figure 
15,  except  that  it  i3  for  the  Hellerman  global  wind  stress 


An  additional  global  data  set  is  the  National 
Meteorological  Center  ( NMC)  analysis  on  2.5°  by  2.5° 
northern  and  southern  hemispheric  grids  where  the  0 
longitude  forms  the  western  boundary.  The  1000-mb  u-  and  v- 
wind  components  from  this  analysis  for  the  period  1  July 
1976  through  29  May  1982  at  twice-daily  intervals  were 
received  from  the  National  Center  for  Atmospheric  Research 
( NCAR)  in  the  NMC  packed  data  format  that  is  described  in 
Office  Note  84  (NMC,  1979).  The  two  hemispheric  grids  were 
merged  into  one  145  by  73.  2.5°  by  2.5°  global  grid  and 
shifted  so  that  1 80W  forms  the  western  boundary.  As  was 
done  for  the  previously  discussed  data  sets,  the  xx  and  ty 
wind  stress  components  were  computed  for  each  analysis  time. 
Interpolation  over  time  (not  more  than  three  days)  filled  in 
the  missing  times  for  which  the  winds  were  not  available. 
After  interpolating,  the  mean  wind  stress  components  were 
computed  for  each  month  of  each  available  year. 

NORAPS 


In  addition  to  the  NOGAPS  model,  FNOC  has  a  Naval 
Operational  Regional  Atmospheric  Prediction  System  (NORAPS) 
for  production  of  numerical  forecasts  over  limited  areas  of 
the  earth.  The  northern  hemisphere  is  one  such  area  for 
which  this  system  is  implemented  on  a  63  by  63  polar 
stereographic  grid  with  the  North  Pole  at  the  center  of  the 
grid  and  a  grid  distance  of  381.00  km  at  60N.  The 
Mediterranean  Sea  63  by  63  polar  stereographic  grid  is  a 
subset  of  this  grid  that  has  been  zoomed  four  times  with  no 
rotation.  It  has  a  grid  distance  of  95.25  km  at  60N.  It  is 
from  the  Mediterranean  grid  that  only  those  u-  and  v-wind 
components  in  the  immediate  vicinity  of  the  western 
Mediterranean  Sea  are  retrieved  from  FNOC  on  a  twice-daily 


basis  and  stored  at  NORDA.  Before  computing  the  x  and  x 

x  y 

wind  stress  components,  the  u-  and  v-wind  components  were 
translated  from  the  polar  stereographic  grid  to  a  latitude- 
longitude  grid  and  interpolated  to  the  109  by  93,  0.1°  by 
0.05°  western  Mediterranean  Sea  Model  staggered  grid,  i.e., 
the  southwest  corner  of  the  u-wind  component  grid  is  at 
35. 2N,  354. 55E  and  that  for  the  v-wind  component  is  at 
35.225N,  354. 5E.  The  wind  stress  components  and  their 
monthly  means  were  then  computed  for  September  1984,  the 
first  complete  month  of  data  received  through  August  1985, 
the  latest  complete  month  available.  Figure  17  shows  the 
wind  stress  vectors  and  Figure  18  the  wind  stress  curl  over 
the  NORDA  western  Mediterranean  Model  region  for  October 
1984  and  January,  April  and  July  1985. 


VI 


STATISTICAL  SURFACE  TO  SUBSURFACE  DATA  CONVERSION 


Hurlburt  (1985)  discusses  the  capabilities  of  a 
numerical  model  of  ocean  dynamics  to  forecast  subsurface 
pressure  information  from  altimetrically  derived  sea  surface 
height  data.  This  project,  intended  to  complement 

Hurlburt's  work,  explores  the  capabilities  of  a  statistical 
model  to  nowcast  subsurface  pressure  information  from  data 
obtained  by  a  hypothetically  perfectly  accurate  model 
altimeter.  Using  height  data  derived  solely  by  a  numerical 
ocean  model,  this  study  employs  a  series  of  regression 
models  to  relate  the  surface  information  to  subsurface 
information. 

The  testbed  region  was  the  Gulf  of  Mexico, 
approximated  by  a  1600  km  by  900  km  rectangular  box  with  20 
km  by  18.75  km  grid  resolution.  The  model  used  was  the  same 
hydrodynamic,  two  layer,  finite  depth,  primitive  equation 
ocean  circulation  model  used  in  Hurlburt's  work.  Selected 
model  simulations  provided  the  "true"  data  for  the 
regressions,  the  predictions,  and  the  error  calculations. 
Experiments  were  run  using  data  from  simulations  with 
barotropic  and/or  baroclinic  instability,  with  a  short 
and/or  a  long  major  time  scale,  and  with  flat  or  large 
amplitude  bottom  topography.  The  two  model  runs  chosen  for 
detailed  study  were: 

El  -  a  nearly  periodic,  barotrop ically  unstable,  run 
with  large  amplitude  bottom  topography,  and 
E2  -  a  highly  baroclinically  unstable  run  with  flat 
bottom  topography. 

Both  simulations  generate  two  pressure  fields,  Pi 
and  P2,  upper  (surface  water)  and  lower  (deep  water) 
pressure,  respectively.  Furthermore,  although  El  and  E2 
have  vastly  different  dynamics,  each  is  characteristic  of  a 
major  dynamical  regime  present  in  the  world's  oceans  (but 
not  necessarily  of  the  actual  Gulf  of  Mexico). 


The  first  series  of  experiments  used  pointwise 
linear  regressions  of  every  PI  point  from  the  model  to  the 
P2  point  directly  beneath,  an  approach  similar  to  that  of 
DeMey  and  Robinson  (personal  communication) ,  who  assumed 
spatial  homogeneity  and  solved  one  regression  equation  for 
the  entire  POLYMODE  region.  Our  results  were  also  similar, 
about  70  percent  spatial  NRMS  error  for  El  and  80  percent 
spatial  NRMS  error  in  E2  for  predictions  made  on  independent 
data  sets.  These  were  clearly  insufficient  for  successful 
deep  water  prediction. 

The  next  series  consisted  of  a  pointwise 
multilinear  regression  of  a  set  of  synoptic  PI  data  to  a  P2 
point.  For  computational  and  statistical  stability  the 
regression  was  done  with  only  the  statistically  meaningful 
principal  components  of  the  PI  set  chosen.  Correlations 
were  used  to  identify  ideal  choices  of  predictors.  The 
tradeoff  is  between  having  sufficiently  independent 
information  in  the  PI  field,  which  often  means  a  greater 
distancing  of  predictors,  and  having  information 
sufficiently  correlated  with  the  P2  field,  which  often 
implies  predictors  crowded  around  the  center  of  the 
prediction;  however,  reasonable  but  not  optimal  choices  of 
predictors  rarely  harmed  the  regressions'  accuracy 
significantly.  Often  a  5  point  star  of  the  PI  point 
directly  above  the  P2  point  being  predicted  and  points 
directly  north,  east,  south,  and  west  of  it,  with  a  radius 
of  40  to  60  km,  was  sufficient. 

This  technique  predicted  El  very  well  on  an 
independent  data  set,  with  spatial  NRMS  error  of  about  0.2 
to  0.3.  Figures  19  and  20,  respectively,  show  "real" 
(model  derived)  lower  layer  pressure,  and  predicted 
(statistically  derived)  lower  layer  pressure  fields,  at  a 
representative  time  in  this  experiment.  As  in  all  the 
figures,  north  is  at  the  top.  All  major  features  are 
accurately  depicted  both  in  strength  and  in  location,  as  is 


evidenced  by  the  low  spatial  NRMS  error  for  the  grid.  The 
multilinear  principal  components  regression  (MPCR)  also  gave 
good  but  increasingly  worse  results  on  sister  experiments 
with  greater  baroclinic  instability  and  less  periodicity. 
In  El  the  eddy  viscosity,  A,  is  300  m2/sec.  A  similar  run 
with  A  -  150  gave  predictions  with  NRMS  of  0.3  to  0.5,  and 
with  A  -  100,  NRMS  is  0.4  to  0.65.  Eddy  mean  energetics  for 
all  experiments  listed  in  this  report  can  be  found  in 
Hurlburt  [1985;  Experiments  T3A  (El),  T3B,  T3C,  and  T1 
( E2 ) ]  . 

On  E2,  errors  were  still  unacceptably  high,  with 
NRMS  of  0.65  to  0.75.  It  was  noticed,  as  might  be  expected, 
that  predictions  were  acceptable  in  mostly  barotropic 
subregions  with  little  eddy  activity,  but  were  very  poor 
otherwise,  often  misplacing  or  totally  misrecognizing  deep 
water  eddies.  Furthermore,  one  dominant  structure  in  this 
region  of  poor  skill  was  the  major  eddy  in  the  upper  layer 
and  an  associated  deep  water  counter-rotating  vortex  pair. 

Figures  21  and  22  show  the  upper  layer  pressure 
field  (PI)  and  the  lower  layer  pressure  field  (P2), 
respectively,  for  a  day  near  the  end  of  the  fifth  year  of 
experiment  E2.  In  Figure  21  an  eddy  is  about  to  be  shed 
from  the  loop  current,  as  a  previously  shed  eddy  begins  to 
die  to  the  west.  Note  the  counter  rotating  vortex  pair  in 
the  lower  layer,  and  the  great  variability  in  the  lower 
layer  pressure  field.  Figure  23  contours  a  prediction  using 
a  MPCR  of  the  PI  field  to  the  P2  field.  Most  of  the  error 
comes  from  severe  inaccuracy  in  the  western  and  southwestern 
portion  of  the  basin,  the  region  of  high  baroclinic 
instability  and  of  frequent  eddy  activity.  In  the 
barotropic  northern  and  northeastern  region  of  the  Gulf  the 
results  are  fairly  accurate. 

The  above  considerations  motivated  the  use  of  a 
limited  Lagrangian  coordinate  MPCR.  A  square  box  the  length 
of  the  major  Loop  Current  eddy  (about  300  km),  and  centered 


on  it  was  allowed  to  move  with  the  eddy  through  a  complete 
eddy  cycle.  The  regression  was  performed  relative  to  the 
box,  selecting  predictors  in  a  five-point  star  as  previously 
stated.  In  E2  on  the  dependent  data  set  the  NRMS  in  the 
region  of  the  eddy  is  0.4  to  0.6,  implying  a  significant 
increase  in  skill  over  the  Eulerian  coordinate  formulations. 
Figures  24  and  25,  respectively,  show  the  "real"  and 
predicted  lower  layer  pressures  in  the  box  described  above 
at  the  same  time  in  E2  as  is  shown  in  Figures  21-23.  The 
region  in  Figure  24  corresponds  to  the  box  drawn  in  the 
right  center  of  Figure  21 .  One  can  see  a  considerable 
increase  in  skill  in  the  prediction,  both  in  determining  the 
strength  and  location  of  the  large  anticyclonic  eddy  to  the 
south  and  in  prediction  of  the  corner  of  the  eddy  to  the 
northwest. 

Kindle  (1985)  notes  the  sufficiency  of  three  or 
four  cross  altimeter  tracks  to  accurately  resolve  an 
irregularly  shaped  eddy.  This  corresponds  in  our 
discretization  to  about  80  km  by  65  km  resolution. 
Subsampling  the  entire  grid  to  this  degree  of  resolution 
still  permitted  statistical  predictions  without  a  serious 
loss  of  accuracy;  that  is,  the  sparser  resolution  usually 
gave  acceptable  results  when  the  finer  resolution  gave 
acceptable  results.  This  intimates  that  practical  usage  of 
the  statistical  model  given  alt imetrically  derived  data  is 
feasible . 

Although  the  Eulerian  methods  tried  show  only 
nowcast  skill,  not  forecast  skill,  since  the  decorrelation 
time  scale  of  each  model  was  only  20  days  or  less,  dynamical 
forecasts  using  the  Eulerian  MPCR  computed  initial  P2  field 
have  often  been  successful  (Hurlburt,  1985;  personal 
communication).  Combining  the  statistical  and  dynamical 
approaches,  using  statistically  derived  information  for 
initial  fields  and  forecasting  an  optimal  estimation  of  the 
P2  field  from  the  statistically  and  dynamically  predicted 
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FIGURE  1:  (a)  Instantaneous  view  of  the  interface  deviation  in  a  two  layer 

Gulf  of  Mexico  model  driven  solely  by  inflow  through  the  Yucatan 
Straits.  The  contour  interval  is  25  m  with  solid  contours  repre¬ 
senting  downward  deviations.  (b)  Depth  of  the  22°C  isothermal 
surface,  4- 1 8  August  1 966  from  Leipper  (1370).  The  contour  in¬ 
terval  is  25  m. 


FIGURE  1 


FIGURE  2:  (a)  Interface  deviation  from  the  simulation  at  day  1970  after  an 
eddy  has  separated  from  the  model  Loop  current  and  propagated 
westward,  (b)  Ninety  days  later  the  eddy  at  day  1970  has  devel¬ 
oped  into  a  counter  rotating  vortex  pair  in  the  western  Gulf. 

The  cyclonic  vortex  is  to  the  north  and  the  anticyclonic  to  the 
south.  The  contour  interval  is  25  m. 


FIGURE  U:  Winter  (a)  vector  plots  of  the  seasonal  climatological  wind 

stresses  from  ship  observations  used  to  drive  the  model,  based 
on  Elliot  (1979).  Maximum  wind  stresses  shown  are  about  1.5 
dynes  in  the  winter.  (b)  Corresponding  steady  state  interface 
deviations  from  the  wind  driven  model  on  a  0.2  degree  grid. 

The  contour  interval  is  12.5  m. 
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FIGURE  5:  Spring  (a)  vector  plots  of  the  seasonal  climatological  wind 

stresses  from  ship  observations  used  to  drive  the  model,  based 
on  Elliot  (1979).  Maximum  wind  stresses  shown  are  about  1.5 
dynes  in  the  winter.  (b)  Corresponding  steady  state  interface 
deviations  from  the  wind  driven  model  on  a  0.2  degree  grid. 

The  contour  interval  is  12.5  m. 


FIGURE  6:  Summer  (a)  vector  plots  of  the  seasonal  climatological  wind 

stresses  from  ship  observations  used  to  drive  the  model,  based 
on  Elliot  (1979).  Maximum  wind  stresses  shown  are  about  1.5 
dynes  in  the  winter.  (b)  Corresponding  steady  state  interface 
deviations  from  the  wind  driven  model  on  a  0.2  degree  grid. 

The  contour  interval  is  12.5  m. 


FIGURE  7-  Fall  (a)  vector  plots  of  the  seasonal  climatological  wind 

stresses  from  ship  observations  used  to  drive  the  model,  based 
on  Elliot  (1979)-  Maximum  wind  stresses  shown  are  about  1.5 
dynes  in  the  winter.  (b)  Corresponding  steady  state  interface 
deviations  from  the  wind  driven  model  on  a  0.2  degree  grid. 

The  contour  interval  is  12.5  m. 


FIGURE  8:  Instantaneous  view  of  the  interface  deviation 
every  90  days,  from  day  90  of  model  year  9  to  day 
0  of  model  year  10,  for  Experiment  28  (left)  and 
Experiment  31*  (right).  Experiment  3^  is  identi¬ 
cal  to  28  except  for  the  addition  of  wind  fo:  c- 
inR.  The  contour  interval  is  25  m. 


FIGURE  9 
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FIGURE  10:  Interface  deviation  mean  and  variability,  for  the 
Gulf  of  Mexico  from  ocean  model.  (a)  Experiment 


31,  wind  forcing  only;  (b)  Experiment  28,  port 


FIGURE  11:  Bottom  topography  and  coastline  geometry  for 
North  Atlantic  model  on  0.25  by  0.2  degree  grid. 
The  contour  interval  is  500  m  and  the  shallowest 
depth  is  1000  m.  The  topography  is  from 
SKNBAPS,  but  has  been  smoothed  twice  to  remove 
small  scale  noise. 


M0DEL  LABEL  =  511 

V.  S.  N0RTH  ATLANTIC  ISN00TH=2) 
DX,DY  =  0.250,0.200  (DEG) 


FIGURE  12:  ( a)  Instantaneous  view  of  the  interface  devotion 

in  a  two-layer  simulation  of  the  North  Atlantic 
driven  from  rest  to  statistical  equilibrium 
solely  by  monthly  climatological  winds.  The 
contour  interval  is  25  m,  with  solid  contours 
representing  downward  deviations. 
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FIGURE  13: 


Bottom  topography  and  coastline  geometry  r°e 
riiif  of  Mexico  model  on  0.1  degree  gri  . 

^  °r  interval  is  500  »  and  the  shallowest 

°d°epth  Is  500  ro.  This  is  the  «w  topography 

field  from  SYNBAPS,  it  is  not  suitable 

direct  use  in  an  ocean  model. 
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FIGURE  14:  Bottom  topography  and  coastline  geometry  for 
Gulf  of  Mexico  model  on  0.1  degree  grid.  The 
contour  interval  is  500  m  and  the  shallowest 
depth  is  500  m.  The  topography  is  based  on  that 
in  Figure  13,  but  is  now  in  a  form  suitable  for 
use  in  the  ocean  model. 


FIGURE  15:  (a)NOGAPS  mean  wind  stress  curl  for  January  1983 

( b )  NOGAPS  mean  wind  stress  curl  for  April  1983 

( c)  NOGAPS  mean  wind  stress  curl  for  July  1983 

(d) NOGAPS  mean  wind  stress  curl  for  October  1983 


FIGURE  1 6 :  (a)Hellerman  climatological  wind  stress  curl  for 
January. 

(b) Hellerman  climatological  wind  stress  curl  for 
April. 

(c) Hellerman  climatological  wind  stress  curl  for 
July. 

(d) Hellerman  climatological  wind  stress  curl  for 
October. 
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FIGURE  17:  (a)NORAPS  mean  wind  stress  for  October  1934 

(b) NORAPS  mean  wind  stress  for  January  1985 

(c) NORAPS  mean  wind  stress  for  April  1985. 

(d) NORAPS  mean  wind  stress  for  July  1985. 
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FIGURE  17(c) 
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FIGURE  17(d) 


FIGURE  18:  (a)NORAPS  mean  wind 

( b )  NORAPS  mean  wind 

( c )  NORAPS  mean  wind 

( d )  NORAPS  mean  wind 


stress  curl  for  October  1984 
stress  curl  for  January  1985 
stress  curl  for  April  1985. 
stress  curl  for  July  1985. 
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FIGURE  18(d) 


FIGURE  20:  The  MPCR  predicted  lower  layer  pressure  in 
experiment  El  at  day  23^0.  The  contour  interval 
is  0.08  m*/s*,  and  labels  are  scaled  by  1000. 
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FIGURE  21 : 


The  "true"  upper  layer  pressure  on  day  1710, 
near  the  end  of  the  fifth  year  of  experiment  E2. 
The  square  box  circumscribes  a  300  km  by  300  km 
region  centered  on  the  major  Loop  Current  eddy. 
The  contour  interval  is  0.5  m*/s2. 
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FIGURE  24:  "Real”  lower  layer  pressure  in  the  region  boxed 
in  Figure  21  surrounding  the  center  of  the  major 
Loop  Current  eddy  on  day  1710  of  experiment  E2. 
The  contour  interval  is  0.1  m2/s2,  and  labels 

are  scaled  by  100. 


FIGURE  25:  MPCR  predicted  lower  layer  pressure  for  the 
region  shown  in  Figure  24,  also  at  day  1710.  The 
contour  interval  is  0.1  m2/s2,  and  labels  are 

scaled  by  100. 


